Structure and melting behavior of classical bilayer crystals of dipoles 
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We study the structure and melting of a classical bilayer system of dipoles, in a setup where 
the dipoles are oriented perpendicular to the planes of the layers and the density of dipoles is 
the same in each layer. Due to the anisotropic character of the dipole-dipole interactions, we find 
that the ground-state configuration is given by two hexagonal crystals positioned on top of each 
other, independent of the interlayer spacing and dipolar density. For large interlayer distances these 
crystals are independent, while in the opposite limit of small interlayer distances the system behaves 
as a two-dimensional crystal of paired dipoles. Within the harmonic approximation for the phonon 
excitations, the melting temperature of these crystalline configurations displays a non-monotonic 
dependence on the interlayer distance, which is associated with a re-entrant melting behavior in the 
form of solid-liquid- solid-liquid transitions at fixed temperature. 
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PACS numbers: 68. 65. Ac Multilayers, 61.50.-f Crystalline state, 52.27.Lw Dusty or complex plasmas; plasma 
crystals, 82.70.Dd Colloids, 64.70.Dv Solid-liquid transitions 
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I. INTRODUCTION 

The realization of a degenerate dipolar gas of 
52 Cr atoms 1,2 and the experimental progress in 
the realization of cold molecular ensembles 3 have 
spurred interest in the properties of particles with 
large dipole moments 4 ' 5,6 ' 7,8,9,10 ' 11,12 in atomic 
and molecular setups 13 . In particular, the strong 
anisotropic dipole-dipole interactions induced in 
ground-state polar molecules by external elec- 
tric fields hold promises for applications ranging 
from the realization of novel strongly-correlated 

p nageg 14, 15, 16, 17, 18, 19, 20, 21, 22. 23, 24, 25, 26. 27, 28, 29, 30, 31, 32, 33, 34 

to quantum simulations 35,36,37 and quantum 
computing 33,38 ' 39,40,41,42 . 

In Refs. 12,29 it is shown that a collisionally stable 
two-dimensional (2D) setup where particles interact via 
purely repulsive effective dipole-dipole interactions can 
be realized by polarizing the molecules using an exter- 
nal electric field, and thus inducing strong dipole-dipole 
interactions, and by confining the motion of the parti- 
cles to a 2D geometry, e.g. by trapping them into a 
single well of a deep optical lattice directed parallel to 
the electric field. The low-energy phase of an ensemble 
of interacting bosonic polar molecules is then a super- 
fluid or a self-assembled crystal for comparatively weak 
and strong interactions, respectively, where the strength 
of the interactions can be modified by varying the inten- 
sity of the polarizing electric field. The crystal is a two- 
dimensional hexagonal lattice structure with quantum 
dynamics given by longitudinal and transverse acous- 
tic phonons. Contrary to familiar Wigner crystals in- 
duced by strong Coulomb interactions 43 , dipolar crys- 
tals emerge at large densities, where dipole-dipole inter- 
actions dominate over the kinetic energy of the particles. 




o » a 



FIG. 1: (Color online) (a) Sketch of the system setup de- 
picting the dipolar bi-layer system: The dipoles are oriented 
along the z direction, d = de z , and confined to two (x,y) 
planes separated by an interlayer distance I along the z di- 
rection. The intralayer dipole-dipole repulsion gives rise to 
a crystalline structure in each layer with basis vectors &i, 
here illustrated by two hexagonal structures with an in-planc 
displacement c. (b) Various possible equilibrium geometries 
for the bilayer crystal, shown as a projection onto the (x, y) 
plane. From left to right: a matching hexagonal (MH) config- 
uration, a one-component hexagonal (OCH) configuration, a 
zigzag rectangular configuration with aspect ratio 02/01 (ZR), 
and a zigzag square (ZS) configuration. The position of the 
dipoles and their nearest-neighbor links in the upper (lower) 
layer are shown by solid (dashed) circles and lines, respec- 
tively. A detailed summary of the lattice parameters is given 
in Table I. 
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This scenario for the realization of 2D crystals can be 
implemented using closed-shell polar molecules as, e.g., 
SrO or RbCs. 

While the scenario above is realized by populating a 
single well of the confining optical lattice, in general it 
will be possible to populate more than one well of the 
lattice. It is thus natural to consider the phases of bi- 
and multilayer configurations of classical and quantum 
dipoles. As a first step, in this work we focus on the 
crystalline structures of a bilayer system of classical 
dipoles. The analogous quantum problem will be the 
subject of a separate study. 



state configuration and we find that the dependence of 
the sound-velocities on the inter-layer separation I is 
non-monotonic, due to the anisotropic character of the 
dipole-dipole interaction. For the excited-states we as- 
sess the stability of the configurations studied in Sec. II 
under phonon fluctuations and find that several of the 
excited states are meta-stable. In Section IV we study 
the thermal melting of the bilayer-crystal. We derive 
the classical melting temperature T m via a modified two- 
dimensional Lindemann criterion, and find that the tem- 
perature T m behaves non-monotonically with increasing 
inter-layer distance. 



In this paper we consider a system of 2N dipoles 
confined into two parallel two-dimensional planes along 
the (x, y)-directions, separated by an interlayer spacing I 
along z. Each plane has the same number of dipoles, N, 
with their dipole-moment d = de z aligned perpendicular 
to the plane, see Fig. 1(a). The interactions between two 
dipoles separated by r are the dipole-dipole interactions 
V(r) = d 2 [l - 3z 2 /r 2 ]/r 3 , where r = |r| > and 
z = r • e z . For our setup this results in the intra-layer 
interactions being always repulsive, which gives rise 
to a hexagonal crystalline structure in each layer. In 
contrast, the inter-layer interactions are attractive 
(repulsive) for dipoles separated by a distance r < y3l 
(r > V3Z), see Fig. 1(a). As detailed in Sees. II-IV, this 
anisotropy in the inter-layer interactions determines that 
the groundstate configuration of the system is a bilayer 
crystal comprised of two hexagonal crystals stacked on 
top of each other. In addition, it gives rise to a non- 
monotonic behavior of the dynamic properties of this 
bilayer crystal (such as the phonon sound velocities and 
the melting temperature) when increasing the interlayer 
distance /, while keeping the intra-layer densities n fixed. 
In particular, within the harmonic approximation for 
the phonon excitations we find that for certain fixed 
temperatures the system displays a reentrant melting be- 
havior in the form of solid-liquid- solid-liquid transitions 
as a function of the dimensionless inter-layer distance 
£ = ly/ri. The picture emerging from these studies is 
one where for large interlayer separations £ 3> 1 the 
two hexagonal crystals of the bilayer structure behave 
as independent crystals, while for vanishing inter-layer 
separations ( < 1 they behave as a single 2D crystal 
of particles with double mass and double dipole-moment. 

The paper is organized as follows: In Section II we 
introduce our model and derive expressions for the po- 
tential energy in the form of a rapidly convergent sum 
obtained via the Ewald summation method. By com- 
paring the energies of a number of possible crystal ge- 
ometries, we determine the structure and energy of the 
ground-state and of several excited-state configurations. 
In Section III we analyze the dynamic properties corre- 
sponding to these crystalline structures within the har- 
monic approximation for the phonon excitations. We 
discuss the phonon excitation branches of the ground- 



II. GROUND AND EXCITED STATES 

We consider a system of classical dipoles confined into 
two planes along the x — y direction, separated by a dis- 
tance I along z. We focus on the situation where the 
dipole-moments di of the particles are aligned perpen- 
dicular to the planes, i.e. along z, and where one has the 
same density n of dipoles in each layer, see Fig. 1(a). The 
interactions between two dipoles is given by their dipole- 
dipole interactions V(r) = d 2 [l — 3z 2 /r 2 ]/r 3 , where 
r = |r| > is their distance and z — r • e z their interlayer 
distance. We note that the dipole-dipole interactions are 
long-range, i.e. decaying like 1/r 3 and anisotropic in 
space. At zero-temperature the system is in a crystalline 
configuration, where the particles in each layer form a 2D 
crystal, and the relative position of the particles is cor- 
related by the long-range dipole-dipole interaction. We 
denote the 2D position of the dipoles within the upper 
(lower) layer by R+j (R_j), which we parameterize by 



R± ,j = iiai + J2a 2 ± c/2, 



(1) 



where the integers j = {31,32) label the j th particle in 
each layer, ai and a 2 are the basis vectors of the peri- 
odic structure with density n and c is a two-dimensional 
vector accounting for a relative in-plane displacement of 
the two structures, see Fig. 1(a). 
The interactions are given by 
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where the first (second) term describes the intra-layer 
(inter-layer) interactions. Since the intra-layer interac- 
tions do not depend on the inter-layer separation i, it is 
convenient to split the energy per dipole, E, into its intra- 
and inter-layer part as E — V/2N = E + Ej. Making 
use of the translational invariance of the infinite system 
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the two contributions read 

( / 2 



2 ^ |R,-|3 ; 



d 2 (|R,- 



E _ f » yi^j t- ^| - 21 ) 
' 2y (|R,+c| 2 +P) 5/2 



(3a) 
(3b) 



where Rj = R CT j — Ro-,0 = ji a i +J2&2 denote the relative 
(2D) positions of the dipoles in each layer. Given the 
slow convergence of the sum in real space involved in 
Eq. (3), we use the Ewald summation method to obtain 
an expression for E involving rapidly convergent sums. 
The explicit derivation is given in Appendix A. We here 
provide only the derived expressions for Eq and Ej , which 
read 
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with |Rj| = (|Rj + c| 2 + Z 2 ) 1 / 2 . In Eq. (4) n denotes the 
intra-layer density and Gj are the 2D reciprocal vectors 
of a single layer, which are parametrized by Gj = jih\ + 
j2^>2 m terms of the primitive translation vectors of the 
reciprocal lattice bi and b2. The quantity a > is an 
(arbitrary) inverse length, for which a convenient choice 
is the inverse of the mean particle separation, i.e., a = 
l/r = 0m. 

In order to determine the ground-state configuration, 
in the following we calculate the energy of a number of 
possible crystal configurations. 

Motivated by the fact that the ground-state for a 
single layer is given by a hexagonal lattice structure, 
we first consider a hexagonal structure in each layer 
with the two structures displaced by c, see Fig. 1(a). 
We find that for an arbitrary inter-layer separation l 7 
the minimal energy is attained for c = (modulo the 
lattice constant). This is also what one would intuitively 
expect, since the closest dipoles in different layers 
tend to attract each other, which leads to a locking of 
the relative positions of the two layers. The attained 
"matching" hexagonal (MH) structure is illustrated in 



TABLE I: Lattice parameters of four considered configura- 
tions: the matching hexagonal (MH), the one-component 
hexagonal (OCH), the zigzag rectangular (ZR) and a zigzag 
square (ZS). From left to right: ai and a.2 are the primitive 
vectors, c is the inter-lattice displacement vector, bi and b2 
are the primitive translation vectors of the reciprocal lattice 
and n is the density of each layer. We introduced the vector 
notation (x, y) = xe x + ye y for the two-dimensional compo- 
nents; the lattice constant is a = |ai| and a-zjax denotes the 
aspect ration for ZR configuration. 

lattice ai/a &ija 2c h\a/2-K b2a/27r no 2 



MH (1,0) a,% 



(L?!) (O'Tf) 75 

OCH (1,0) (0,^3) (1,V3) (1,0) (0,^) ^ 

ZS (1,0) (0,1) (1,1) (1,0) (0,1) 1 

ZR (1,0) (0,£) (1,£) (1,0) (0,£) £ 



Fig. 1(b) and its basis and reciprocal vectors are listed 
in Table I, together with three other (mcta-stable) 
configurations detailed below. In Figure 2 we plot 
the energy per particle i5^ MH ^ attained for the MH 
configuration at a fixed density n as a function of the 
dimensionless inter-layer separation ^=Z/(y / 7rro), given 
by the inter-layer distance in units of the mean particle 
separation (up to a constant l/y^r), as a solid line. We 
notice that for large inter-layer separations, £ ^> 1, the 
ground-state energy approaches the value corresponding 
to a single hexagonal (SH) layer in its ground-state, 
i.e. £( MH )(£ -» 00) = £( SH ) « 4.443d 2 n 3 / 2 . In 
the opposite limit, £ <C 1, the energy is domi- 
nated by the large attraction between the two layers 

£(MH)(£ < 1) ~ (_1/£3 + 2 x 4 .443) rf 2 n 3/2 ; which di _ 

verges as cx — l/l 3 for I — > 0, due to the strong attraction 
of the closest dipoles in different layers. We remark 
that the latter corresponds to an unphysical regime 
for typical molecular systems 12 . In fact, at these short 
distances the required depth of the optical trapping 
potential for realizing a 2D layer becomes unreasonably 
large. Moreover one expects short-range interactions 
in the form of, e.g., Van-dcr-Waals interactions as well 
as the core repulsion, to dominate the interactions at 
these small distances 44 . These latter effects have been 
neglected in our model, c.f. Eq. (2). However, we notice 
that here the physically relevant quantity is £, and the 
limit £ — > can be approached by decreasing the density 
of dipoles n while keeping I fixed. 

In order to confirm that the MH configuration is the 
ground-state, we compare its energy to those of a number 
of other (intuitively motivated) configurations 45 . 

As a first alternative configuration we consider a "one 
component" hexagonal (OCH) structure, which is ob- 
tained by removing every second dipole in each layer in 
a staggered way and rescaling the relative density. The 
OCH configuration is illustrated in Fig. 1(b), and its ba- 
sis and reciprocal vectors arc summarized in Table I. The 
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FIG. 2: (Color online) The energy per dipole E for a fixed 
density n as a function of £ for three different lattice configu- 
rations: MH (solid line), OCH (dash-dotted line) and ZS (dot- 
ted line) . The inset is a blow-up of the excited state energies in 
the region of small interlattice separation, < £ < 0.3, show- 
ing a crossing of the OCH- and ZS-energies around £ as 0.12. 
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FIG. 3: (Color online) The energy per particle E as a function 
of £ for the OCH (dash-dotted line), the optimal ZR (dashed 
line) and ZS configuration (dotted line), showing the transi- 
tions: OCH — > ZR — > ZS with increasing £. Inset: The aspect 
ratio a-z/a-i for the ZR-configuration as a function of £. 



attained energy 

^(OCH) ig plotted in Fig 2 function 
of £ as a dashed-dotted line, and we notice that it exceeds 
that of the MH configuration. The label OCH has been 
chosen for this configuration, as it resembles a hexagonal 
lattice for a single component, when looked from above. 
Accordingly, in the limit of vanishing inter-lattice sepa- 
ration the energy of the OCH configuration tends to that 
of a single hexagonal layer (with a double density), i.e. 

£(OCH) ( £ = 0) „ 4 44 3 d 2 (2n) 3/2 M 12 57 g ^3/2 

As a second alternative configuration we consider a 
"zig-zag" square (ZS) structure, where dipoles in each 
layer form a square lattice, and the two layers are shifted 
with respect to each other as illustrated in Fig. 1. The ba- 
sis and reciprocal vectors for the ZS structure are given in 
Table I. The energy for the ZS-structure E^- zs ^ is shown 
in Fig. 2 as a dashed line. The figure shows that for 
large £ the energy of the ZS configuration exceeds the 
MH energy but it is smaller than the OCH energy, while 
at small £ the energies of the ZS and OCH configurations 
become comparable. The inset of Fig. 2 is a blow-up of 
the excited state energies in the region < £ < 0.3. From 
the inset, we observe that the two energies, £(° CH ) and 
E( zs \ actually cross at £ = £ ~ 0.12 and for £ < £ the 
OCH structure is energetically favored over the ZS one. 

The excited state configurations OCH and ZS can be 
interpolated smoothly by "stretching" the lattice, which 
in general leads to a "zig-zag" rectangular (ZR) struc- 
ture with a variable aspect ratio 02/01 in each layer. 
The ZR structure is illustrated in Fig. 1(b) and its ba- 
sis and reciprocal vectors are listed in Table I. In order 
to find the optimal ZR configuration, we minimize the 
energy as a function of the free parameter 02/01- The 
obtained aspect ratio 02/01 for the optimal ZR config- 
uration in shown in the inset of Fig. 3, and we notice 



that the resulting structure coincides with the OCH and 
ZS configurations at £ = and £ « 0.17, respectively. 
The energy obtained for the optimal ZR configuration, 
E( ZR \ is plotted as a dashed line in Fig. 3 in the region 
0-05 < £ < 0.17 , along with the ones obtained for the 
OCH (dash-dotted line) and ZS structures (dotted line). 
The figure shows that E (z ^ equals £( OCH ) and E (zs ^ at 
£ = and £ as 0.17, respectively, while it is lower in the 
parameter region in between. 

The results above are consistent with the MH config- 
uration being the ground-state of the system for all £, 
due to the attraction between particles in the two layers 
(separated by r < This in contrast to the situ- 

ation occurring for Coulomb bilayer systems, where the 
repulsion between particles in the different layers leads 
to a change of the ground-state configuration depending 
on the ratio of the inter-layer separation to the mean 
intra-layer spacing 45 . However, we notice that analogous 
(smooth) transitions between different configurations oc- 
cur here between excited-state structures, which will be 
shown below to be metastable. These structures may in 
principle be prepared using properly-designed in-plane 
optical lattice potentials, as argued below. 



III. DYNAMICAL PROPERTIES 

In this section, we study the phonon spectra for those 
equilibrium lattice configurations, which we found in the 
previous section. Thereby we make use of the harmonic 
approximation of the lattice excitations, and determine 
the stability (meta-stability) of the obtained configura- 
tions for the ground (excited states) under small fluctu- 
ations. 



5 



In order to obtain the excitation spectra for the various 
lattice configurations, we consider the dynamical matrix 
Mc(q), whose eigenvalues are the square of the phonon 
frequencies 46 . We notice that the bilayer configurations 
above correspond to a (single 2D) Bravais lattice with a 
unit cell given by two molecules, one in the upper and 
one in the lower layer. Thus Md (q) is here a 4 x 4 matrix 



are given by 



M D (q) 



/ Ml x + 


M xy 


Mf_ 


M xy _ \ 


M x + \ 


M vv 


M x + V _ 


M v + V _ 


M xx_ 


M xy _ 


M* x _ 


M xy _ 


V M+L 


M v + V _ 


M x _ v _ 


M y _ y _ ) 



(5) 



where the superscript r = x (r = y) refer to a displace- 
ment along x (y) and the subscript a = + (a = — ) refer 
to the components in the upper (lower) layer. The matrix 
elements in Eq. (5) read 



M rv 

(7(7 



Ml 



1 

m 
1 



[DY + (Q) - DY + ( q ) + DY_(0)] , (6a) 



-#;-(q)], 



(6b) 



where m is the mass of the dipoles. The quantities 
D™ a (q) are defined as 



d\ u M) = ^e-"i-( Ro '+- R -)a T a^ 0+J(T (r = o),(7) 

3 



with Vo+j a (r) the two-body interaction potential be- 
tween the dipole at position in layer + and the dipole 
j in layer a. At q = the quantities D^(q) for a = + 
(a - = — ) correspond to the intra-layer (inter-layer) force 
constants. 

Using the Ewald summation method [see Eq. (A9) and 
Eq. (Af7) in Appendix A], we can rewrite the sum in 
Eq. (7) into the following rapidly convergent forms 



D 7 + (<1) = E(q + Gj)r(q + G^T (^±31,0 

3 

+ — S TV + V lim d T d v Qi + r|) , 



3*0 



(8a) 



D 7M) = £(q + G,) T (q + G,), e ^- c T(i^±3i ; a i 

3 

+ e- i *< R t+ c ) lim drdM2(\Rj + r\) , 



(8b) 



where 5 TV is the Kronecker-delta, and where we used the 
notation (q + Gj) T = e r • (q + Gj) for the component in 
the direction r of the two-dimensional (reciprocal) vec- 
tors. The functions T(x,y), f2i(a;) and ^(x) in Eq. (8) 



fii(x) 



-in 



-x -y 



+ J2 ( ±2 ) uxe ±2xy eric (x±y), 

(9a) 
(9b) 
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2a e~ a x 



erfc(ax) 2a e a 



3l 2 



erfc(aa;) 2a (3 + 2a 2 x 2 ) e 



Z^pK x 4 



(9c) 



For the configurations of Fig. f , the complex Hcrmitian 
matrix Mo(q) of Eq. (7) can be transformed into a real 
and symmetric matrix. This is achieved by first applying 
the unitary transformation Mo(q) = UMd^US with 
U a 4 x 4 matrix defined as 



(10) 



with I2 the 2x2 identity matrix. This transformation 
brings the dynamical matrix into the symmetric form 




M D (q) 



/ M ++ + Im[M^ 
= ^ Rc[M + _] 

with M + „ ee 



M 



Rc[M+_] 
Im[M+_] 



M xx 

M x + l 



++ 



MH 



M. 



mi 



(11a) 
(lib) 



That the matrix Mrj(q) is real now stems from the fact 
that Im[M-| ] vanishes for a lattice with inversion sym- 
metry 45 ' 47 , which is the case for all lattice configurations 
considered in this work. 

For each quasi-momentum q, diagonalizing Mrj(q) 
provides the square of the phonon frequencies, ui a (q) 2 
(with 1 < a < 4), which correspond to the four dis- 
tinct phonon modes of the bilayer system. Within the 
harmonic approximation, the stability of the various lat- 
tice configurations is linked to the sign of the eigenvalues 
of Mrj(q). That is, a lattice configuration is stable if 
all four eigenvalues of Md (q) are positive for all quasi- 
momenta q, i.e. u> a (q) 2 > 0, while it is unstable if one 
(or more) of the four eigenvalues is negative for a given 
quasi-momentum q, c.f. o; Q (q) 2 < 0. 



A. Ground state configuration 

In this subsection we calculate the phonon spectrum 
for the ground-state configuration MH as a function 
of the inter-layer separation £, using the techniques 
described above. As expected, we find both acoustic and 
optical modes, with the latter related to out-of-phasc 
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FIG. 4: (Color online) Phonon dispersion curves w«(q) for 
the MH lattice configuration in units of ujo = yj d 2 n 5 / 2 /m, 
for £ = 0.1 (solid lines), £ = 0.7 (dashed lines) and £ — > oo 
(dash-dotted lines). The frequencies are presented along the 
high symmetry directions in the Brillouin zone. The high- 
symmetry points r, X and J are depicted in the inset. Note 
that the different axis scaling for the high-frequency regime 
U) 3> 6a;o . 



vibrations of particles in different layers. We show 
that the longitudinal and transverse sound velocities of 
the acoustic modes show a non-monotonic dependence 
on £, which is linked to the anisotropic nature of the 
dipole-dipole interaction. The picture emerging from 
these studies is one where the bilayer structure behaves 
for vanishing inter-layer separations £ <C 1 as a single 
2D crystal of particles with double mass and double 
dipole-moment, while for £ 1 the two layers behave as 
independent 2D crystals. Within the harmonic approxi- 
mation inherent to the present discussion, the transition 
between these two situations occurs approximately for 
inter-layer distances such that the nearest-neighbor 
inter-layer interaction switches from repulsive (( < 1) 
to attractive (£ 1). 

Figure 4 shows the phonon dispersion relations for the 




+ 0.5 1.0 1.5 2.0 

% 

FIG. 5: (Color online) Sound velocities for the MH lattice 
configuration as a function of the dimensionless layer separa- 
tion £: (a) Longitudinal sound velocity wla and (b) transverse 
sound velocity vta, in units of u)o/y/n. 



MH lattice configuration along the high symmetry direc- 
tions in the Brillouin zone for a three values of £, i.e. 
£ = 0.1,0.7,100 (solid, dashed, dash-dotted lines). The 
frequencies are in units of the characteristic phonon fre- 
quency ojo = V d 2 n 5 / 2 /m. The symmetry points T, X 
and J are depicted in the inset. The figure shows that 
for £ — ► oo there are two phonon branches. Each one 
of these branches is doubly degenerate, corresponding to 
independent longitudinal and transverse acoustic phonon 
modes of the two layers. We find that, as expected, the 
modes exactly match those of a single layer, confirm- 
ing that in the limit £ — > oo the two layers behave as 
independent. For finite £ the dipole-dipole interaction 
couples the two layers, and lattice vibrations in the two 
layers become correlated. Accordingly, the figure shows 
that for £ = 0.7 (dashed lines) and £ = 0.1 (solid lines) 
the phonon modes develop separate acoustic and optical 
branches. Since the MH lattice structure can be repre- 
sented as the repetition of a basis cell comprising two 
particles, one per layer, stacked on top of each other, 
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FIG. 6: (Color online) Optical frequencies uj op (in logarith- 
mic scale in units of ljo) at the V point for the MH lattice 
configuration as a function of £ (solid line), along with the 
approximations ln(w p/wo) ~ 0.716 — 2.502 ln^ for £ < 0.7 
(dashed line) and ln(w op /wo) w 3.919 - 3.441£ for £ > 0.7 
(dotted line). 



the optical modes are easily understood as arising from 
out-of-phase vibrations of the dipolcs in each basis cell. 
The figure shows that the optical frequency of vibration 
increases with decreasing inter-layer distance, and the 
bandwidth of the optical branch tends to shrink (note 
that in the figure the scales for the optical modes for 
£ = 0.7 and £ = 0.1 arc different). In the (unphysi- 
cal) limit £ — > 0, the phonon spectrum reduces to one 
where there is only one, i.e. non-degenerate, longitudinal 
acoustic and one transverse acoustic mode, while the op- 
tical branches are shifted to infinitely large frequencies. 
This observation is consistent with the system behaving 
as a single layer of dipoles with double mass and double 
dipole strength, arranged in a hexagonal configuration 
(see below). 

Figure 4 shows that the acoustic branches for £ = 0.7 
have lower frequencies than the corresponding branches 
for £ = 0.1 and £ — > oo, indicating a non-monotonic 
dependence of the frequencies on £. In order to bet- 
ter investigate this point, in Fig. 5(a) and Fig. 5(b) 
we plot the longitudinal and transverse sound veloci- 
ties v LA = du)L A /dq |g =0 , «ta = duj TA /dq \ q=0 , with 
wla and wta the frequencies of the longitudinal and 
transverse acoustic modes, respectively. We find that 
for £ 3> 1, the longitudinal and transverse sound veloci- 
ties tend to i>la — 2.544wo/v / "' an d ^ta — 0.768wo/\/"-> 
respectively. These values correspond to the sound ve- 
locities of a monolayer hexagonal lattice configuration, 
consistent with the observation above that for £ 3> 1 the 
crystal vibrations in the two layers become independent. 
In the opposite limit £ — > 0, wc find that the sound ve- 
locities are larger than those above exactly by a factor 
y/2. A simple comparison with the characteristic phonon 



frequency ojq = \J d 2 n 5 / 2 /m shows that this v^-factor is 
consistent with the system behaving as a single layer of 
dipoles of mass 2m and dipole strength 2d. 

Figure 5 shows that the dependence of the sound ve- 
locities on £ is non-monotonic. In particular, minima 
for «la and Vta occur at £ = £o « 1 and 0.64, respec- 
tively (see Appendix B). This non-monotonic behavior 
is linked to the anisotropic character of the dipole-dipole 
interaction. In fact, for £ — > oo the two layers behave as 
independent hexagonal crystals. For finite £, the inter- 
layer interactions couple the two layers. This coupling, 
which is responsible for the formation of the optical band, 
splits the degeneracy of the phonon frequencies, lowering 
the acoustic sound velocity, which is associated with a 
softening of the crystal, see Fig. 4. For £ < 1 the optical 
and acoustic branches are well separated, corresponding 
to the formation of a crystal of tightly bound pairs of 
dipoles. For £ — > the sound velocity is a \/2 larger 
than at £ — > oo, as discussed above, and this determines 
the appearance of a minimum in between, and we ob- 
serve that the value £ « 1 roughly corresponds to the 
inter-layer distance at which the interaction of a dipole at 
R+j with next-nearest-neighbor in the opposite layer at 
position R-j' switches from attractive to repulsive. We 
remark that although the crystalline structure is softened 
for finite £, the sound velocity remains always finite. This 
is due to the anisotropic character of the dipole-dipole 
interactions, which ensures that the MH configuration is 
always the stable ground-state configuration. 

Figure 6 shows as a solid line the optical frequencies 
LJop as a function of £ in logarithmic scale at the T point, 
where the two optical frequencies are degenerate. We 
notice that for £ > 0.7 the frequencies decay exponen- 
tially with increasing £ as ln(w op /a;o) ~ 3.919 — 3.441£ 
(c.f. dotted line), while they diverge as a power law 
\n(uj op /uj ) w 0.716 - 2.502 ln(£) for £ < 0.7 (c.f. dashed 
line). This change in behavior is another manifestation 
of the crossover from two independent layers to a single 
crystal of paired dipoles. 



B. Excited-state configurations 

In this subsection wc study the meto-stability of the 
excited-state configurations OCH, ZS and ZR intro- 
duced in Sect. II, by calculating the phonon modes for 
each configuration. For a given interlayer distance £, 
within the applicability of the harmonic approximation, 
regions of meta-stability and instability for the above 
configurations correspond to real and imaginary values 
of the computed phonon frequencies, respectively. By 
analyzing the sound velocities of the phonon excita- 
tions, we find that the instability of the OCH and ZS 
crystalline structures is associated with the vanishing of 
the transverse acoustic branch of the phonon modes in 
the directions TX and TM. We thus derive a stability 
diagram for the low-lying excitations of the system as 
a function of £, see below. This is interesting, since 
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FIG. 7: (Color online) Square of the phonon frequencies a> 2 
in units of uiq = d 2 n 5//2 /m (a) for the OCH configuration at 
£ = (solid lines), £ = 0.2 (dashed lines) and £ — > oo (dash- 
dotted lines) and (b) for the ZS configurations at £ = 0.1 
(solid lines), £ = 0.2 (dashed lines) and £ — + oo (dash-dotted 
lines). The frequencies are shown along the high symmetry 
directions in reciprocal space for each lattice configuration. 
The high-symmetry points F, X and M are depicted in the 
insets. 



FIG. 8: (Color online) (a) Longitudinal sound velocities vla 
and (b) transverse sound velocities vta in units of ujo/y/n 
along the (1,0) (solid lines) and the (1,1) (dotted lines) direc- 
tions for OCH, ZR and ZS lattice configurations. The vertical 
dotted line denotes the boundary between ZR and ZS. The 
gray shading corresponds to the transition region where the 
harmonic approximation for the phonon excitations may be- 
come inadequate. 



in principle these excited-state configurations may be 
realized, e.g. by first trapping cold polar molecules in 
optical lattices with the same geometry as OCH and ZS 
crystals, increasing the dipole-dipole interactions using 
external fields, forming interaction-induced OCH and ZS 
crystals, and finally adiabatically removing the lattice 
potential. 

Panels (a) and (b) in Fig. 7 show the square of 
the phonon frequencies for the OCH and ZS lattice 
configurations, respectively, for a few values of £ and 
in units of uJq = d 2 n 5 / 2 /m. The phonon spectra in 
Fig. 7 are shown along the high symmetry directions 
in reciprocal space, with symmetry points labeled in 
the insets. The figure shows that for certain values of £ 
the square of the phonon frequencies becomes negative, 
lo 2 < 0, signaling an instability of the corresponding 
crystalline structure for the given value of £. In panels 



(a) and (b), regions of meta-stability are present for 
£ = and £ = 0.2, respectively. In fact, we determined 
numerically that the OCH configuration of panel (a) is 
meia-stable for £ < 0.170, while the ZS configuration of 
panel (b) is ?neia-stable in the range 0.166 < £ < 0.247. 
We notice that in these regimes where lo 2 > the system 
in the OCH and ZS configurations is meta-stable, since 
the associated crystalline structures are excited states 
of the system. By analyzing the sound velocities of 
the phonon excitations, we found that the instability 
of the OCH and ZS crystalline structures is associated 
with the vanishing of the transverse acoustic branch of 
the phonon modes in the directions (1,1) and (1,0), as 
detailed below. 

Panels (a) and (b) in Fig. 8 show the sound velocities 
of the transverse acoustic (TA) and longitudinal acous- 
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tic (LA) modes for the OCH, ZS and ZR configurations 
as functions of £ in the range of stability of each con- 
figuration, respectively. In the figure, the continuous 
and dashed lines correspond to the (1,0) and (1,1) di- 
rections, respectively. Panel (a) shows that the TA mode 
for the OCH configuration vanishes in the (1,0) direc- 
tion at £ ~ 0.170 [sec also Fig. 7(a)]. The TA mode 
for the ZS configuration vanishes in the (1,1) and (1,0) 
directions for £ < 0.166 and £ > 0.25, respectively In 
addition, panel (b) shows that the longitudinal modes 
for all three configurations soften with increasing £. The 
transverse and longitudinal sound velocities for the ZR 
configuration interpolate between those of the OCH and 
ZS configurations, see Fig. 8(a) and Fig. 8(b), coinciding 
with those at £ = and £ ~ 0.166, respectively. 



monic approximation for the excitations of the crys- 
tal using a (modified) Lindcmann criterion. The latter 
states that for a given configuration the melting occurs 
when the mean relative displacement between neighbor- 
ing sites becomes of the order of the mean interparticle 
distance ro = 1/ y 7 ™ (see Table. I for the configurations 
of Fig. I) 48 . This reads 



y/([u(R)-u(R + r)] 2 ) 



OCH 



ZS >i 



MH 



0.0 



0.166 0.170 0.247 



FIG. 9: Approximate regions of stability for the MH, ZS, and 
OCH configurations, as a function of £. 

Figure 9 summarizes the regions of stability for the MH 
configuration, which is the ground-state for any £ > 0, 
the OCH and the ZS configurations . We remark that the 
values £ =0.166, 0.170 and 0.247 have been obtained nu- 
merically, in the harmonic approximation for the phonon 
spectrum. Since the harmonic approximation is bound 
to break down around any transition points between var- 
ious configurations, the numerical values 0.166 and 0.170 
should be taken with caution, and simply interpreted as 
indicative of the transition region. 

For completeness, in Fig. 10 we show the optical fre- 
quencies at the T point as function of £ for the OCH, ZR 
and ZS lattice configurations. As expected, we find that 
the frequencies corresponding to the ZR configuration 
interpolate between those of the OCH and ZS configu- 
rations. In particular, they become degenerate around 
£ w 0.170, where the ZS configuration becomes stable. 
This degeneracy is a natural consequence of the fact that 
the aspect ratio aijay of the ZS configuration is 1. We 
notice that the different behaviors of the optical frequen- 
cies may be used to distinguish experimentally the vari- 
ous metastablc structures, when initially prepared in tai- 
lored optical lattice potentials, as discussed at the begin- 
ning of this subsection. 



IV. CLASSICAL MELTING 

In this section we discuss the classical melting tem- 
perature of the bilayer system, as obtained in the har- 



where (|u(R) — u(R + r)| 2 ) is the relative mean square 
displacement, u(R) and u(R + r) are the displacement 
vectors at site R and at its nearest neighbor site R + r, 
() is a thermal average, and 5 m < 1 is a parameter which 
in general has to be determined numerically 



The left-hand side of Eq. (12) is computed in the har- 
monic approximation for the phonons as follows. Each 
particle in the bilayer structure has two different kinds 
of nearest-neighbors: in-plane and out-of-plane. Thus, in 
analogy to the Coulomb case of Ref. 45 we define the in- 
tralayer ( Au ++ ) and interlayer ( Am_| ) correlation func- 
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FIG. 10: (Color online) Optical frequencies uj op in unit of ojq 
at the r point for the OCH, ZR and ZS lattice configurations. 
The vertical dotted line denotes the boundary between ZR 
and ZS. The gray shading corresponds to the transition region 
where the harmonic approximation for the phonon excitations 
may become inadequate. 
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FIG. 11: (Color online) Illustration of the dipole- 
configuration for the MH structure leading to the modified 
Lindemann criterion. 



tions (see Appendix C) 
Au ++ 

2k B T 



j- E E <K(o)-<wi 2 > 

+ f=x,yh=l...S + 



^E^-y{[^(q^-) 2 + P, + (q,.) 2 ] 

[l-cos(q-R+)] j, (13a) 
7T E E H(0)-u-(h)\ 2 ) 



~f=x,y h=l...S- 

k B T x - 



+Ay(q,j) 2 + Pyi^jf -2 p+(ci,j)p-(q,j) 



+Py(lJ)Py(<lJ) cos(q-R h 



(13b) 



Here is the number of nearest-neighbor dipoles in 
layer a = ±, u°(h) is the 7 th component of the displace- 
ment of a particle at position /i in the layer a, p°{<\,j) 
is the 7 th component of the eigenvector of the j th mode 
at point q in the Brillouin zone of the sublattice in layer 
a, and is the relative vector connecting one particle 
to its h nearest-neighbor in the same (t = +) or op- 
posite (r = — ) layers. We notice that the number of 
nearest-neighbors S a and their distance depends on the 
considered lattice configuration (see Fig. 1) and £. In 
particular, Fig. 11 shows that, in the relevant case of the 
groundstate configuration MH (sec Sect. IV A below), the 
number of in-plane nearest-neighbors is 6, while that of 
out-of-plane nearest-neighbors is 7. 

The correlation ([u(R)-u(R + r)] 2 ) of Eq. (12) is now 
computed as 

([u(R)-u(R + r)] 2 ) = Au+++f(l)Au+-, (14) 



where the function f(l) describes the influence of lattice 
vibrations in one layer on the lattice vibrations in the 
opposite layer, and it is defined as 



f(l) 



1 



-3kP 



(1 + ^ 2 ) 5 / 2 {1 + kP) 7 /2 



(15) 



Here, the geometric parameter n can be obtained from 
Table I, and it reads k = (a 2 ) -1 and k = ( | c | 2 ) 1 for the 
MH and zigzag lattice configurations, respectively This 
expression for /(/) is chosen to be proportional to the 
in-plane part of the force between two nearest-neighbor 
dipoles in opposite layers and it satisfies the conditions 



lim/m 



1 



and lim f(l) 

l — >oo 



0, 



where the latter condition is due to the fact that vibra- 
tions in the two layers are independent for infinite inter- 
layer separations. 



A. Melting of the ground state configuration 

In this subsection we determine the classical melting 
temperature T m of the ground-state crystal configuration 
MH as a function of £, using the modified Lindemann 
criterion introduced above. We find a non-monotonic 
dependence of T m on £, which we attribute to the 
anisotropic nature of the dipole-dipolc interactions. 
This is interesting since for certain temperatures it is 
associated with a re-entrant melting behavior in the 
form of solid-liquid- solid-liquid transitions. 

Figure 12 shows the melting temperature T m as a func- 
tion of £, as calculated from the Lindemann criterion 
with S m = 0.23 within the harmonic approximation for 
the phonon modes. The precise value of S m should in 
principle be obtained numerically, e.g. using molecular 
dynamics simulations. In the absence of such a computa- 
tion for a classical bilayer crystal, the value of S m = 0.23 
has been chosen in analogy to the one obtained in Ref. 30 
for the quantum melting transition from a single-layer 
crystal of bosonic dipoles into a superfluid using Dif- 
fusion Monte-Carlo techniques. Using this value of 8 m 
in Eq. (12), we find that for £ ^> 1 the classical melt- 
ing temperature of the bilayer crystal tends to the value 
To « 0.066d 2 /a 3 . By construction, the latter corresponds 
to the classical melting temperature of a single hexagonal 
crystal as computed in the harmonic approximation dis- 
cussed above. We here notice that the obtained value of 
T is of the order of the actual one T num w 0.089d 2 /a 3 for 
a classical single-layer crystal, as obtained numerically by 
molecular dynamics simulations 49 . Since the spirit of the 
Lindemann criterion is that of a qualitative estimate of 
the transition point, in the following we will be content 
with the value To. 

For ^ <C 1 the figure shows that the melting tempera- 
ture tends to T m = 4To. This is consistent with the pic- 
ture of a hexagonal crystal made of paired dipoles with 
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FIG. 12: (Color online) (a) Melting temperature T m as a function of £ for the MH lattice configuration (solid line). The unit 
To ~ 0.066d 2 /a 3 corresponds to the melting temperature of a single-layer crystal, as computed from Eq. (12) with S m = 0.23. 
The dashed line is a guide to the eye. (b) Melting temperature T m as a function of £ for particles interacting via the potential 
Vlj^ of Eq. (16). Here, r is the strength of the attractive part of the potential (see text), with t = 1 corresponding to the 
dipole-dipole interaction of panel (a). The dotted line is a guide to the eye. 



dipolc strength d* = 2d, as discussed in Sect. Ill A. Inter- 
estingly, the figure shows that T m has a non-monotonic 
dependence on £ around £ « 1. In particular, the T m - 
vs-£ curve has a local maximum and a local minimum at 
£ « 1 and £ f=a 0.6, respectively. As noticed in Sect. Ill A, 
this region of ^-values corresponds to the distance at 
which the dipole-dipole interaction between a particle in 
one layer and its nearest-neighbor in the opposite layer 
changes sign from attractive to repulsive (e.g. d$ and d~^ 
in Fig. 1, respectively, with 1 < j < 6). 

In order to check that the non-monotonicity is in fact 
connected with the anisotropic nature of the dipole- 
dipole interaction, in Fig 12(b) we have plotted the melt- 
ing temperature for an artificial system of particles where 
the strength of the attractive part of the dipole-dipole in- 
teraction can be tuned, and thus the particles interact via 
a potential of the form 

Vg* = d 2 ( -L + T Zlll) (16) 

with r a constant, < r < oo, and | r^- 1 the intcrparticle 
distance. We find that the non-monotonic character of 
the curve is enhanced for r > 1, while it tends to disap- 
pear for t < 1 and in particular it vanishes for t < 0.9. 
In the limit r — > of purely repulsive interactions (not 
shown) the system resembles the Coulomb case of Ref. 45 , 
and accordingly we find that the MH lattice ceases to be 
the groundstate configuration. 

The observations above confirm that the reentrant 
(non-monotonic) behavior of T rn as a function of £ is due 
to the attractive character of the dipole-dipole interac- 
tions. This indicates that it is possible to alternate solid 
and liquid phases by changing the intcrlaycr distance or 



the density of the molecules. 



V. CONCLUSION 

In this work we studied the structure, the stability 
and the melting of a classical bilayer system of dipoles, 
polarized perpendicular to the layers, as a function of 
the interlayer distance and the density of dipoles in 
each layer. Using the Ewald summation technique, we 
have computed the ground-state energy and the phononic 
spectrum of a few physically-motivated lattice config- 
urations, finding that the ground state is always the 
matching hexagonal crystal configuration, where two tri- 
angular lattices are stacked on top of each other, as ex- 
pected. Higher-energy configurations have been found to 
be metastable in different regimes of interlayer distances 
and dipole densities. These configurations may be real- 
ized using polar molecules trapped in optical lattices of 
properly chosen geometry. 

The main result of this work is that the anisotropic na- 
ture of the dipole-dipole interaction potential profoundly 
affects the dynamical properties of the bilayer system of 
dipoles, and its melting behavior. In fact, on one hand 
we found that the attractive part of the potential deter- 
mines a non-monotonic dependence on £ of the longitu- 
dinal and transverse sound velocities in the ground-state 
configuration. On the other hand, we have shown that 
the classical melting temperature of the bilayer crystal 
has an interesting reentrant behavior as a function of £. 
This reentrant behavior is due to the anisotropy of the 
dipole-dipole interaction, and it entails that it is possi- 
ble to alternate various crystalline and liquid phases by 
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changing £ at a fixed temperature. 

The present analysis is motivated by the recent devel- 
opments in the physics of cold molecular gases, which 
may provide for a physical realization of these sys- 
tems. In particular, the classical melting of the bilaycr 
crystalline phases may be realized in future in setups 
where cold polar molecules are trapped in adjacent wells 
of a one-dimensional optical lattice, and their dipoles 
arc polarized by a static electric field oriented paral- 
lel to the optical lattice. Under appropriate trapping 
conditions 12 ' 29 , the resulting in-plane dipole-dipole in- 
teractions are purely repulsive, while inter-plane interac- 
tions can be repulsive or attractive. 
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with Rj = Rcrj — R CTj o and |Rj + r| = (|Rj + c + r| : 
l 2 ) 1 ' 2 . Then E and Ei can be obtained from 



E Q = limrf 2 Vo(r,0), 

r— >0 

Ei = limd 2 V/(r,0). 

r— >0 



(A3a) 
(A3b) 



We use the identity based on the integral representation 
of the gamma function 



1 1 

x 2 * ~ f(s) 



/ t s_1 cxp(-x 2 t)dt, 
Jo 



(A4) 



with s = 3/2,r(3/2) = 0F/2 for ip and ip n , s = 
5/2,r(5/2) = 3Vtt/4 for ip I2 , and the 2D Poisson sum- 
mation formula 

^exp(-|p + R 3 | 2 ;-zq- (p + R,)) 

3 

^ t" 1 ^ cx P (iG, ■ P ) exp (-^±Sl!) (A5) 



L 



where Gj = j±h\ + j2^2 (with integers i,j ) is the two- 
dimensional vector in reciprocal lattice, L 2 = 1/n is the 
area per primitive cell. Then tpoi^j^l) can be expressed 
by 



APPENDIX A: RAPIDLY CONVERGENT FORM 

OF Vo AND ipi 

The direct numerical computation of sums over lat- 
tice sites with long-range dipole-dipole interactions is 
in general impractical. Thus in the following we trans- 
form each sum into rapidly convergent forms using Ewald 
method. 45,50 The detailed techniques are shown in the 
following. For the calculation of the rapidly convergent 
form of the energy, at first we define the following two 
functions 



-0o(r,q) 
M r . q) 



e -*q-(Rj+r) 



l R j 



(Ala) 



iq-(Rj+c+r) 



3 



|R j+ r|3 

_3 ^2 giq^Rj +c+r) \ 



where 
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E 
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E 



iq-(Rj+c+r) 
iq-(Rj+c+r) 



|R,- 



(Alb) 

(A2a) 
(A2b) 



exp 



|G i+ q| 2 \,. 2 



dt - 



9 f 00 

Ye-in-n-i' / i^HR.+rl 2 ^ (A6) 



3^0 

where a is a small positive number. After using the in- 
tegration 

/>OC 

/ t 1 / 2 cxp(-\x\ 2 t)dt 

Ja 2 

aexp(— a 2 |x| 2 ) 



2\x\- 



and 



erfc(a|2:| 



t -l/2 exp /_ l±L. )dt 



(A7) 



xlxl^crfc ( — ) }, 



(A8) 



where the expression contains the complementary error 
function erfc(z) = 1 — erf(z) = J* e~* dt, we ob- 
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tain the final form of tpo(r, q) 



In the same way, we transform ipi2 as 



-2\Gj +q| erfc 



exp 
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4a 2 
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Similarly, 

^ = ^E e * (qH 



f 2a\ 3 + 2a 2 |R,+r| 2 „_ a 2| R ,. +r |S 
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(A14) 



-1/2 



where the integrations 



x exp ( ^ Gj 4 * q|: - l 2 t) dt + -±= E e-^^ +c ) 



i 1/2 exp (-|Rj + r| 2 t) dt, (A10) 



by replacing t by w 2 , the first integration of above equa- 
tion can be rewritten as 
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and 



Using the integration 
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and Eq. (A7) we have the form of ijjii 



,i(q+G,-)T JG r c 
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2a 



are used. After simplifying, we obtain the rapid conver- 
gent form of ipi(r, q) as 
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FIG. 13: (Color online) The quantities F£*, Fyy 1 , Fll and 
Fyy as functions of £. 
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Therefore, the rapid convergent forms of Eq and Ei are 
expressed as Eqs. (4a) and (4b). 



APPENDIX B: INTERPRETATION OF THE 
MINIMA OF SOUND VELOCITIES FOR MH 
CONFIGURATION 

For understanding the minimum of the LA and TA 
modes of the sound velocities of MH configuration, we 
define the following functions 
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Where Ej is the interaction potential between a dipole 
at origin and the dipoles in the opposite layer, u x r y -\ is 
the x(y) component of the displacements of the origin 
dipole around its equilibrium position, r x t y -\ is the x(y) 
component of the distance between the origin dipole and 
the dipoles in the opposite layer. 

The first two plots in Fig. 13 describe and Fyy 
as functions of £, which indicate that the direction of 
vibration is along the direction of propagation. We find 



that the minimum is located at the same value of £ in 
the curves of and Fyy as that in the plot of ula [see 
Fig. 5(a)], where £ m j n w 1. The remaining two plots of 
Fig. 13 show Fxx and Fy* vs £, which express that the 
direction of vibration is perpendicular to the direction of 
propagation, the value of £ for the minimum in this two 
curves is the same as that in vta [see Fig. 5(b)], where 
£min = 0.64. Due to the intrinsic property of interlaycr 
dipole-dipole interaction, the four functions above have 
minima at £ m j n « 1 and £ m i n = 0.64 respectively. 



APPENDIX C: CORRELATION FUNCTIONS 
FOR BILAYER SYSTEM 

From the Fourier transformation we know that the 7 th 
component of the displacement vectors of a dipole at the 
origin position in layer + and at the h th nearest neighbor 
position in layer a are 



1,3 

^^ C ^(q,jX(q,j)exp(zq.R^),(Cl) 



<(0) 



<(h) 



where p^(q, j) is the 7 component of the eigenvector of 
j th mode at q point in the first Brillouin zone of the 
sublattice in layer a = ±. c^(q, j) is the probability 
parameter of pZ(q,j), RJJ is the relative position of the 
h th nearest neighbor dipole in layer a. Making use of the 
relation (c^(q,j)c^(q',j')> = (fc B 7 1 /mw 2 (q, j)) 8 qtq >6 jd ,, 
we can obtain the relative mean square displacements 
between the two considered nearest neighbors 

(|«+(0)-<(/0| a ) 



2k B T 
Nm 



I? ^by { w (q ' j)2] (1 - cos q • R ^ ) } ' 



(C2a) 



(\u+(0)-u-(h)\ 2 ) 



- 2 [p+(q,i)p~(q, j)] cos(q-R^) \, 



(C2b) 



where u~(0) and itj 



\h) are the 7 component of the 
original and the h th nearest neighbor dipoles in the +(— ) 
layer, m is the mass of the dipoles, ks is the Boltzmann 
constant, w(q, j) is the phonon frequency of j th mode at 
q point in the first Brillouin zone. After the summation 
over the 7 th components and over the nearest neighbor 
sites, we can finally obtain the expressions Eqs. (13a) 
and (13b). 
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